install.packages("rgdal", repos = "https://packagemanager.posit.co/cran/2023-10-13")1 Overview
Network constrained Spatial Point Patterns Analysis (NetSPAA) is a collection of spatial point patterns analysis methods special developed for analysing spatial point event occurs on or alongside network. The spatial point event can be locations of traffic accident or childcare centre for example. The network, on the other hand can be a road network or river network.
In this hands-on exercise, we will gain hands-on experience on using appropriate functions of spNetwork package:
to derive network kernel density estimation (NKDE), and
to perform network G-function and k-function analysis
2 The Data
In this study, we will analyse the spatial distribution of childcare centre in Punggol planning area. For the purpose of this study, two geospatial data sets will be used. They are:
| Type | Name | Details |
|---|---|---|
| Geospatial | Punggol_St |
|
| Geospatial | Punggol_CC |
|
3 Installing and launching the R packages
In this hands-on exercise, four R packages will be used, they are:
| Package | Description |
|---|---|
| sf | Provides functions to manage, processing, and manipulate Simple Features, a formal geospatial data standard that specifies a storage and access model of spatial geometries such as points, lines, and polygons. |
| spNetwork | Provides functions to perform Spatial Point Patterns Analysis such as kernel density estimation (KDE) and K-function on network. It also can be used to build spatial matrices (‘listw’ objects like in spdep package) to conduct any kind of traditional spatial analysis with spatial weights based on reticular distances. |
| tidyverse | A collection of functions for performing data science task such as importing, tidying, wrangling data and visualising data. |
| tmap | For thematic mapping; provides functions for plotting cartographic quality static point patterns maps or interactive maps by using leaflet API |
To install and launch the four R packages.
pacman::p_load(sf, spNetwork, tmap, tidyverse)4 Import Data and Preparation
4.1 Import
The code chunk below uses st_read() of sf package to important Punggol_St and Punggol_CC geospatial data sets into RStudio as sf data frames.
network <- st_read(dsn="data/geospatial",
layer="Punggol_St")Reading layer `Punggol_St' from data source
`C:\kytjy\ISSS626-GAA\Hands-on_Ex\Hands-on_Ex03\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 2642 features and 2 fields
Geometry type: LINESTRING
Dimension: XY
Bounding box: xmin: 34038.56 ymin: 40941.11 xmax: 38882.85 ymax: 44801.27
Projected CRS: SVY21 / Singapore TM
childcare <- st_read(dsn="data/geospatial",
layer="Punggol_CC")Reading layer `Punggol_CC' from data source
`C:\kytjy\ISSS626-GAA\Hands-on_Ex\Hands-on_Ex03\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 61 features and 1 field
Geometry type: POINT
Dimension: XYZ
Bounding box: xmin: 34423.98 ymin: 41503.6 xmax: 37619.47 ymax: 44685.77
z_range: zmin: 0 zmax: 0
Projected CRS: SVY21 / Singapore TM
4.2 Checking Contents
networkSimple feature collection with 2642 features and 2 fields
Geometry type: LINESTRING
Dimension: XY
Bounding box: xmin: 34038.56 ymin: 40941.11 xmax: 38882.85 ymax: 44801.27
Projected CRS: SVY21 / Singapore TM
First 10 features:
LINK_ID ST_NAME geometry
1 116130894 PUNGGOL RD LINESTRING (36546.89 44574....
2 116130897 PONGGOL TWENTY-FOURTH AVE LINESTRING (36546.89 44574....
3 116130901 PONGGOL SEVENTEENTH AVE LINESTRING (36012.73 44154....
4 116130902 PONGGOL SEVENTEENTH AVE LINESTRING (36062.81 44197....
5 116130907 PUNGGOL CENTRAL LINESTRING (36131.85 42755....
6 116130908 PUNGGOL RD LINESTRING (36112.93 42752....
7 116130909 PUNGGOL CENTRAL LINESTRING (36127.4 42744.5...
8 116130910 PUNGGOL FLD LINESTRING (35994.98 42428....
9 116130911 PUNGGOL FLD LINESTRING (35984.97 42407....
10 116130912 EDGEFIELD PLNS LINESTRING (36200.87 42219....
str(network)Classes 'sf' and 'data.frame': 2642 obs. of 3 variables:
$ LINK_ID : num 1.16e+08 1.16e+08 1.16e+08 1.16e+08 1.16e+08 ...
$ ST_NAME : chr "PUNGGOL RD" "PONGGOL TWENTY-FOURTH AVE" "PONGGOL SEVENTEENTH AVE" "PONGGOL SEVENTEENTH AVE" ...
$ geometry:sfc_LINESTRING of length 2642; first list element: 'XY' num [1:2, 1:2] 36547 36559 44575 44614
- attr(*, "sf_column")= chr "geometry"
- attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA
..- attr(*, "names")= chr [1:2] "LINK_ID" "ST_NAME"
st_crs(network)Coordinate Reference System:
User input: SVY21 / Singapore TM
wkt:
PROJCRS["SVY21 / Singapore TM",
BASEGEOGCRS["SVY21",
DATUM["SVY21",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4757]],
CONVERSION["Singapore Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (N)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (E)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Cadastre, engineering survey, topographic mapping."],
AREA["Singapore - onshore and offshore."],
BBOX[1.13,103.59,1.47,104.07]],
ID["EPSG",3414]]
childcareSimple feature collection with 61 features and 1 field
Geometry type: POINT
Dimension: XYZ
Bounding box: xmin: 34423.98 ymin: 41503.6 xmax: 37619.47 ymax: 44685.77
z_range: zmin: 0 zmax: 0
Projected CRS: SVY21 / Singapore TM
First 10 features:
Name geometry
1 kml_10 POINT Z (36173.81 42550.33 0)
2 kml_99 POINT Z (36479.56 42405.21 0)
3 kml_100 POINT Z (36618.72 41989.13 0)
4 kml_101 POINT Z (36285.37 42261.42 0)
5 kml_122 POINT Z (35414.54 42625.1 0)
6 kml_161 POINT Z (36545.16 42580.09 0)
7 kml_172 POINT Z (35289.44 44083.57 0)
8 kml_188 POINT Z (36520.56 42844.74 0)
9 kml_205 POINT Z (36924.01 41503.6 0)
10 kml_222 POINT Z (37141.76 42326.36 0)
str(childcare)Classes 'sf' and 'data.frame': 61 obs. of 2 variables:
$ Name : chr "kml_10" "kml_99" "kml_100" "kml_101" ...
$ geometry:sfc_POINT of length 61; first list element: 'XYZ' num 36174 42550 0
- attr(*, "sf_column")= chr "geometry"
- attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA
..- attr(*, "names")= chr "Name"
st_crs(childcare)Coordinate Reference System:
User input: SVY21 / Singapore TM
wkt:
PROJCRS["SVY21 / Singapore TM",
BASEGEOGCRS["SVY21",
DATUM["SVY21",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4757]],
CONVERSION["Singapore Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (N)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (E)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Cadastre, engineering survey, topographic mapping."],
AREA["Singapore - onshore and offshore."],
BBOX[1.13,103.59,1.47,104.07]],
ID["EPSG",3414]]
5 Visualising Geospatial data
par(bg = '#E4D5C9')
plot(st_geometry(network))
plot(childcare,
add=T,
col='#800200',
pch = 19)
par(bg = '#E4D5C9')
plot(network)
plot(childcare,
add=T,
col='#800200',
pch = 19)
For achieving a visually appealing and interactive representation of geospatial data, the tmap package’s mapping function can be employed.
tmap_mode('view')
tm_shape(childcare) +
tm_dots(col = "#800200") +
tm_shape(network) +
tm_lines()6 Network Constrained KDE (NetKDE) Analysis provided in spNetwork
6.1 Preparing the lixels objects
Prior to computing NetKDE, it is necessary to partition the SpatialLines object into lixels with a specified minimum distance. This operation can be accomplished using the lixelize_lines() function from the spNetwork package.
lixels <- lixelize_lines(lines = network, #<< SpatialLinesDataFrame
lx_length = 700, #<< Length of a lixel
mindist = 350) #<< Minimum length of a lixellixels_750 <- lixelize_lines(lines = network, #<< SpatialLinesDataFrame
lx_length = 750, #<< Length of a lixel
mindist = 375) #<< Minimum length of a lixelLixelize_lines Function
- Dimensions for Lixels Objects:
Set the length of a lixel (lx_length) to 700m.
Set the minimum length of a lixel (mindist) to 350m.
After cut, if the final lixel is shorter than the minimum distance, it will be added to the previous lixel.
Segments that are already shorter than the minimum length are not modified.
If the minimum length is NULL, then mindist = maxdist/10.
Additional Information about
Lixelize_linesFunction:Lixelize_linesis used to cut a SpatialLines object into lixels with a specified minimal distance.The function also supports multicore processing through
lixelize_lines.mc().
- Post-cut Considerations:
- After cutting, if the length of the final lixel is shorter than the minimum distance, it is added to the previous lixel.
- If the minimum distance is NULL, then mindist is set to maxdist/10.
- Segments that are already shorter than the minimum distance are not modified.
6.2 Generate line centre points using lines_center() of spNetwork
lines_center()of spNetwork is used to generate a SpatialPointsDataFrame with line center points.Points are located at center of line based on the length of the line.
samples <- lines_center(lixels)samples_750 <- lines_center(lixels_750)6.3 Performing NetKDE
To compute the NetKDE:
densities <- nkde(network,
events = childcare,
w = rep(1, nrow(childcare)),
samples = samples,
kernel_name = "quartic", #<<
bw = 300,
div= "bw",
method = "simple", #<< Can be simple, discontinuous, continuous
digits = 1,
tol = 1,
grid_shape = c(1,1),
max_depth = 8,
#agg = 5, #<< Aggregate events within a 5m radius (faster calculation)
sparse = TRUE,
verbose = FALSE)densities_750 <- nkde(network,
events = childcare,
w = rep(1, nrow(childcare)),
samples = samples_750,
kernel_name = "quartic", #<<
bw = 300,
div= "bw",
method = "simple", #<< Can be simple, discontinuous, continuous
digits = 1,
tol = 1,
grid_shape = c(1,1),
max_depth = 8,
#agg = 5, #<< Aggregate events within a 5m radius (faster calculation)
sparse = TRUE,
verbose = FALSE)Kernel Method and Arguments:
The code chunk reveals the use of the quartic kernel (kernel_name argument).
spNetwork supports various kernel methods, including triangle, gaussian, scaled gaussian, tricube, cosine, triweight, epanechnikov, or uniform.
Calculation Methods for NKDE:
The method argument indicates the use of the “simple” method for calculating NetKDE.
spNetwork offers three methods for NKDE:
simple: Distances between events and sampling points are replaced by network distances. The kernel formula is adjusted to calculate density over a linear unit instead of an areal unit.
discontinuous: Proposed by Okabe et al. (2008), this method “divides” the mass density of an event at intersections of lixels.
continuous: An alternative version proposed by Okabe et al. (2008) adjusts the density before intersections, making the function continuous.
User Guide Reference:
- The user guide of the spNetwork package provides a comprehensive discussion of
nkde(). It is recommended to read the guide to understand various parameters for calibrating the NetKDE model.
Additional Notes on Arguments:
The chosen kernel method is quartic, and the decision is explained.
spNetwork supports alternative kernel methods such as triangle, gaussian, scaled gaussian, tricube, cosine, triweight, epanechnikov, or uniform.
The selected method for NKDE calculation is “simple,” and the reasons for its use are explained.
Other supported methods include “discontinuous” and “continuous,” each with specific characteristics described in the code chunk.
6.4 Visualising NetKDE
6.4.1 Insert computed density values (i.e. densities) into samples and lixels objects as density field
samples$density <- densities
lixels$density <- densitiessamples$density_750 <- densities_750
lixels$density_750 <- densities_7506.4.2 Rescale density values from number of events per m to number of events per km
As the svy21 projection system is in meters, the resulting density values are very small (e.g., 0.0000005). The code below employed to rescale the density values from the number of events per meter to the number of events per km.
# rescaling to help the mapping
samples$density <- samples$density*1000
lixels$density <- lixels$density*1000samples$density_750 <- samples$density_750*1000
lixels$density_750 <- lixels$density_750*1000`
6.4.3 Using tmap package to plot map after rescaling
tmap packages can be used to prepare interactive and high cartographic quality map visualisation.
tmap_mode('view')
tm_shape(lixels)+
tm_lines(col="density")+
tm_shape(childcare)+
tm_dots()tm_shape(lixels)+
tm_lines(col="density_750")+
tm_shape(childcare)+
tm_dots()Interpretation
- Road segments with relatively higher density of childcare centres (darker color)
- Road segments with relatively lower density of childcare centres (lighter color)
7 Network Constrained G- and K-Function Analysis
Objective: Conducting CSR test using the
kfunctions()function from the spNetwork package.Null Hypothesis (\(H_0\)): The observed spatial point events (i.e., distribution of childcare centres) exhibit a uniform distribution over a street network in Punggol Planning Area.
CSR Test Assumption:
The CSR test relies on the assumption of a binomial point process.
Assumption implies that childcare centres are randomly and independently distributed over the street network.
Interpretation of Results:
If the null hypothesis is rejected:
Inference: The distribution of childcare centres shows spatial interactions and dependence.
Resultant Patterns: Nonrandom patterns may be observed.
CSR Test Execution:
- Execution involves utilizing the kfunctions() function from the spNetwork package.
kfun_childcare <- kfunctions(network,
childcare,
start = 0,
end = 1000,
step = 50,
width = 50,
nsim = 50,
resolution = 50,
verbose = FALSE,
conf_int = 0.05)Insights
Arguments Used:
Ten arguments are employed in the code chunk, namely:
- lines: A SpatialLinesDataFrame with sampling points.
- points: A SpatialPointsDataFrame representing points on the network.
- start: Start value for evaluating the k and g functions.
- end: Last value for evaluating the k and g functions.
- step: Jump between two evaluations of the k and g functions.
- width: Width of each donut for the g-function.
- nsim: Number of Monte Carlo simulations (50 simulations in the example).
- resolution: Resolution for simulating random points on the network.
- conf_int: Width of the confidence interval (default = 0.05).
- For additional arguments, refer to the user guide of the spNetwork package.
Output of kfunctions():
The function outputs a list with the following components:
- plotkA: ggplot2 object representing k-function values.
- plotgA: ggplot2 object representing g-function values.
- valuesA: DataFrame with values used to build the plots.
We can visualise the ggplot2 object of k-function by using the code chunk below.
kfun_childcare$plotk +
labs(title ="K-Function") +
theme(panel.background = element_rect(fill = "#E4D5C9"),
plot.background = element_rect(fill = "#E4D5C9"),
panel.grid.major = element_line(colour = "#efe7df", linetype = 1, linewidth = 0.5),
panel.grid.minor = element_line(colour = "#efe7df", linetype = 1, linewidth= 0.5),
plot.title = element_text(face = "bold", size = 12, hjust = 0.5),
)
Observations from the Graph:
- The blue line indicates the empirical network K-function for childcare centers in Punggol.
- A gray envelope represents results from 50 simulations spanning the 2.5% to 97.5% interval.
Inference:
- Blue line values between 125m-687.5m fall below the gray envelope.
Conclusion:
- Implies that childcare centers in Punggol exhibit a regular pattern within the 125m-687.5m distance range.
Kfunctions values:
8 Reference
Kam, T. S. Network Constrained Spatial Point Patterns. R for Geospatial Data Science and Analytics. https://r4gdsa.netlify.app/chap07.html